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A new technique is proposed for the recovery of optical phase from intensity information. The method is 
based on the decomposition of the transport-of-intensity equation into a series of Zernike polynomials. An 
explicit matrix formula is derived, expressing the Zernike coefficients of the phase as functions of the Zernike 
coefficients of the wave-front curvature inside the aperture and the Fourier coefficients of the wave-front 
boundary slopes. Analytical expressions are given, as well as a numerical example of the corresponding phase 
retrieval matrix. This work lays the basis for an effective algorithm for fast and accurate phase retrieval. 



1. INTRODUCTION 

The problem of phase retrieval from intensity informa- 
tion is of relevance in many areas of science. 1,2 In this 
paper we consider the phase retrieval technique based on 
the transport-of-intensity equation (TIE) first proposed by 
Teague. 3,4 We are motivated to investigate this problem 
as the TIE has been studied in the context of adaptive op- 
tics for astronomy 5,6 and the imaging of phase objects in 
microscopy. 7 Related problems have also benefited from 
the sorts of approach described here, such as determina- 
tion of the aberrations of the eye in ophthalmology, 8 cor- 
recting optics for x-ray sources, 9 and the investigation of 
aberrations in large optical telescopes, 10 

In this paper we are concerned not with the nature 
of a specific application but rather with the detailed 
mathematical basis for the solution of the TIE. We 
concentrate here on solutions in terms of the Zernike 
polynomials, 11 as these are the natural starting point 
for the discussion of the diffraction theory of aberrations 
and so are widely used in adaptive optics and other opti- 
cal studies. 12 " 19 Although the details of our discussion 
revolve around these polynomials, it is no doubt also pos- 
sible to develop closely related approaches with the use of 
other orthonormal polynomials, which may be useful in 
specific ^applications. 

The TIE is a partial differential equation that directly 
relates the phase distribution in the planes orthogonal 
to the optical axis to the rate of change of the wave- 
front intensity of the beam. The equation forms the ba- 
sis of the now widely used wave-front curvature sensing 
technique proposed by Roddier and Roddier. 5,6,10 Here 
we suggest a new approach based on Zernike decompo- 
sition of the TIE that we hope will allow phase to be 
sensed with improved speed and resolution. We propose 
to expand each function involved in the TIE into a se- 
ries of Zernike polynomials, thus reducing the boundary- 
value problem for the TIE to a system of linear algebraic 
equations. Such an approach is effective in the case of 
circular apertures with uniform illumination, to which we 
confine our present study. The analysis of the structure 
of the resulting algebraic system significantly clarifies 
the contribution of each particular phase aberration to 
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the propagation of the wave-front intensity distribution 
in the beam. Furthermore, the reduction of the TIE to 
a system of linear equations allows simple and efficient 
methods for its solution^ Some facts about the algebraic 
properties of the TIE with respect to Zernike polynomials 
have been reported earlier. 10 " 19 We give rigorous proofs 
for the relevant results, derive some new ones, and bring 
them together to present a complete picture of the struc- 
ture of the TIE with respect to individual Zernike aber- 
rations of the phase. 

Our study is based oh the calculation of the Laplacian 
of Zernike polynomials. These results can be considered 
as a logical extension of the work by Noll, 13 who calculated 
the first-order partial derivatives of Zernike polynomials. 
We prove that the only Zernike polynomials with zero 
Laplacian are the diagonal ones, i.e., those with radial 
degree equal to the azimuthal frequency. We also prove 
that the Laplacian of any nondiagonal Zernike polynomial 
of radial degree N can be represented as a linear combi- 
nation of Zernike polynomials each with radial degree not 
exceeding N - 2. Using these results, we derive an ex- 
plicit operator formula for the phase solution to the TIE 
as a function of Zernike coefficients of the wave-front cur- 
vature inside the aperture and the Fourier coefficients of 
the wave-front slopes at the boundary. It turns out that, 
on account of the stability of the solution to the Neumann 
problem for the Poisson equation, this phase solution is 
insensitive to small errors in the data of the wave-front 
curvature and boundary slopes. Furthermore, the coef- 
ficients of the inverse operator do not depend on the ex- 
perimental data. Therefore we believe that our approach 
lays the basis for a very efficient algorithm for phase re- 
trieval by the TIE method. • 

In Section 2 we review the basics of the phase recon- 
struction with the TIE. In Section 3 we develop our new 
approach, and we discuss the results in Section 4. 

2. TRANSPORT-OF-INTENSITY EQUATION 
AND PHASE RETRIEVAL 

The underlying idea of phase retrieval with the use of 
the TIE is that in the paraxial approximation the evolu- 
tion of the intensity distribution in the direction of beam 
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propagation is defined mainly by the distribution of the 
phase in the planes orthogonal to that direction. There- 
fore these , phase distributions can be recovered if the 
intensity change from one such plane to another is mea- 
sured. In this section we recall the basics of the intensity 
transport method, as suggested by Teague 3,4 and Roddier 
and Roddier. 5 - 610 

Let us consider the scalar monochromatic electromag- 
netic wave with complex amplitude 



exp(ikz)u(r) = / 1/2 (r)exp[t*2 + i<j)(r)], 



(1) 



where r = (jc, y, z). In the paraxial (Fresnel) approxi- 
mation with the optical axis parallel to z, the complex 
amplitude u(r) satisfies the paraxial equation 



(2ikd z + A)u(x, y, z) = 0, 



(2) 



where k is the wave number, d z = d/dz> and A = V 2 = 
$x 2 + d y 2 is the two-dimensional Laplacian. If we substi- 
tute Eq. (1) into Eq. (2) and separate real and imaginary 
parts, we obtain the following pair of equations [provided 
that /(r) * 0]: 



2kd 2 <f> = -\V<j>\ 2 + Z>(/), 
kd z I = -V/ • - /A<£, 



(3) 
(4) 



where V = (d z , d y ) is the gradient operator in a plane and 
D(I) = / -l/2 A(/ 1/2 ) is the diffraction term. Equation (4) 
is the TIE. It can be used for the reconstruction of the 
phase in some area CI of a plane (x, y), z = constant, 
if the distributions of intensity and its z derivative are 
known there. 

In this work we will consider only circular domains fl , 
where R is the radius and T is the boundary of Ct . It is 
convenient to introduce the polar coordinates (r, 6) in the 
plane of interest, z = 0. We also restrict our study to the 
case of uniform intensity distributions in ft : 



/(r, 6) = I 0 H(R - 
H(t) = 



r), 

1 

0 



Iq = constant , 
f >0 
t <0 



(5) 



The assumption of uniform intensity is widely accepted 
in adaptive optics. 5,6 * 20 Substituting Eqs. (5) into Eq, (4), 
we obtain 

-H(R - r)A<Mr, S) + S(R - r)d r <f>(R, 6) 

= kl 0 - l d z l(r,6) t (6) 

where 6(r) is the Dirac delta function and d r <f> is the 
phase derivative along the radial direction. Equation (6) 
implies that the z derivative of intensity must also contain 
a delta-function term at the boundary: 



kh~ l dj(r t 0) = f(r t 0) + S(R - r)ijt(6), 



(7) 



where the function f is smooth up to the boundary and 
tp is a smooth function on the boundary T. Comparing 
Eqs. (6) and (7), we find that 



-A0-/* 

inside the circular domain H and that 



dr<t> = <A 



(8) 



(9) 



(10) 



on the boundary T. 



Thus, in the case of uniform intensity, the phase can 
be obtained as a solution to the Neumann boundary-value 
problem (9) for the Poisson equation (8). This approach 
was developed by Roddier 5,6 and has an important advan- 
tage over the original suggestion of Teague, 3,4 who consid- 
ered the Dirichlet boundary conditions <f> = tp on T. The 
advantage of boundary conditions (9) is in the fact that, 
unlike the phase itself, the phase normal derivative at 
the boundary can be found from intensity measurements 
(7) at this boundary. Hence direct measurements of the 
phase boundary values are not necessary. 

When one studies a boundary problem for a partial 
differential equation, it is always necessary to address 
three major questions, namely, those concerning exis- 
tence, uniqueness, and stability of solutions. Because 
the Neumann problem [Eqs. (8) and (9)] is a classi- 
cal object of mathematical physics, its properties are 
well known. 

It is proved in the theory of partial differential equa- 
tions that a solution to the problem of Eqs. (8) and (9) 
exists if and only if the following condition holds 21 : 

//a * )rdrd ^ + / *W)*d0 =0. 
When we substitute for/" and i^, Eq. (10) becomes 

k \ dj(r, 0)rdrdd = -I 0 R / d r <f>{R, 6)d6 , 

Jo Jo Jo 

(11) 

which is just an expression of conservation of energy; loss 
of intensity in a region arises through energy flow across 
the boundary of the region. Equation (10) may be used 
to check the consistency of acquired intensity data. 

The mathematical theory also states 21 that the solution 
of Eqs. (8) and (9) is unique up to a constant, i.e., if <f> is 
a solution, then <f> + C is also a solution for any constant 
C. This arbitrary additive constant is not essential for 
the phase reconstruction. A nontrivial fact is that in the 
case of uniform intensity (5) and circular domain fl the 
phase reconstructed by Eqs. (8) and (9) is unique (up 
to a constant) even in the class of multivalued phase 
functions. 22 This is important in view of the example 
given by Gori et al., 23 which presents essentially different 
(multivalued) phase functions corresponding to the same 
(nonuniform) three-dimensional intensity distribution in 
a wave field. 

Finally, we would like to recall that the solution 
<t> to Eqs. (8) and (9) is stable with respect to small 
errors in f or ^, as a result of the boundedness of the 
inverse operator. 21,24 Namely, if <f> and <f>' are the solu- 
tions to Eqs. (8) and (9) with the right-hand-side functions 
if, iff) and (/\ respectively^with d n {f t f) < $\ and 
dri*!* f < &2f where d ft and dr are the appropriate met- 
rics inside H and on the boundary T, respectively, then 
dn(<i>y 4>') < € and e = €(S lt S 2 )— 0 when S x + S 2 — 0. 

Note that, in the approach described above, the 
boundary values of the phase normal derivative ${6) = 
d r <f>(R t $) should be obtained as the coefficient of the 
delta function in Eq. (7). In reality, the intensity change 
near the boundary always has a finite gradient. If the 
intensity is almost uniform in the interior of H and has a 
sharp decrease near the boundary, then at the boundary 
the first term on the right-hand side of Eq. (4) is much 



1934 J. Opt. Soc. Am. A/Vol. 12, No. 9/September 1995 



Gureyev et al. 



larger than the second one, which allows us to write 



= -k 



(12) 



Inside O the first term on the right-hand side of Eq. (4) 
is much smaller than the second one, which gives us the 
expression for f : 



f{r,d) = kl 0 - x d z l(r, 0), 



r <R. 



(13) 



Thus we have a well-defined boundary-value problem 
[Eqs. (8) and (9)] with the right-hand-side functions /"and 
(A obtainable from the measurements of optical intensity 
in two closely spaced planes (we need to measure inten- 
sity in two planes in order to calculate the derivative dj). 
We now proceed with the solution of Eqs. (8) and (9) by 
the method of orthogonal expansions. 

3. TRANSPORT-OF-INTENSITY EQUATION 
AND ORTHOGONAL POLYNOMIALS 

A. Expansion of the Transport-of-Intensity 
Equation into Orthogonal Polynomials 

In this subsection we briefly outline the scheme of the 
orthogonal expansion method of solution of boundary- 
value problems for partial differential equations. In the 
following subsections we will apply it to the TIE using 
Zernike polynomials, though it is possible to implement 
this method with any complete set of orthogonal func- 
tions. In particular, it may be interesting to consider 
the eigenfunctions of the Laplacian in the circle. Our 
choice of Zernike polynomials is motivated by their favor- 
able properties with respect to the description of phase 
aberrations. 11-14 

Let {Zj} be a complete set of linearly independent func- 
tions in domain fl, so that we can expand the phase 
<f> and the wave-front curvature /"into a series over Zy. 
<t> — Y.<f>jZj, f — Z fjZj- Substituting these into the 
Poisson equation (8) and using the linearity of the Laplace 
operator, we obtain the system of linear algebraic equa- 
tions Y.^ij<f>j = fi t where the matrix elements A tJ are 
the coefficients of the decomposition of -&Zj over Z[\ 
(-A)Zj = XA^Z,-. If the system {Zj\ is orthonormal 
with respect to some scalar product, {Z it Zj) = 5,7, where 
5,-y is the Kronecker symbol, then A ,7 = (-AZ;, Z{). If 
the matrix [A/,] is singular, i.e., if it maps some of 
the Zj or their linear combinations into zero, then the 
phase 4> should be expressed as a sum of two components 
(projections), <f> = 4> <0) + <f> {l) , where tf> {0) is the projection 
of <p onto the subspace Ker(- A) spanned by all linear com- 
binations of Zj mapped by the Laplacian into zero, and 
= - 0(0) i s the projection onto the complementary 
subspace [Ker(-A)] 1 . 

By definition {-&)<f> i0) = 0; hence <£ (0) cannot be found 
from the above system, because XA ( 7<£, = X A;^)- 1 ' = fi 
does not. depend on <£ (0) . However, if the boundary prob- 
lem is well posed, 21 <t> {0) can be uniquely found from the 
boundary conditions. The matrix [A,,] is always non- 
singular (invertible) on the subspace [Ker(-A)] 1 ; hence 
the component <f> {l) can be obtained by the inversion of this 
matrix: 4>j li = Y.^ji~ l fi- Thus the phase can be re- 
constructed as a sum of two components, <j> = 0 (O) + 
0 a \ with the component <f> ili obtained from the wave- 



front curvature and the component <f> {0i determined from 
boundary conditions. 

In practice, we must deal with truncated series, which 
is equivalent to considering finite-dimensional subspaces 
spanned by subsets of the whole system {Zj}. In what 
follows, we define such natural subsets and implement 
the method of orthogonal expansions with the system of 
Zernike polynomials. 

B. Zernike Decomposition of the Laplacian 

In this subsection we derive the decomposition of the 
Laplacian in the spaces Z# of Zernike polynomials with 
radial degree not exceeding some integer N. The main 
aim is to find the kernel (the polynomials mapped into 
zero) and the image (the functions into which Zernike 
polynomials are mapped) of the Laplacian in Z N . Such 
an analysis is a necessary preliminary step for the phase 
reconstruction by the Zernike expansion of the TIE, which 
we describe in Subsection 3.C. Examples of the sub- 
spaces that we introduce in this subsection can be found 
in Table 1. 

We recall the definition of Zernike polynomials using a 
notation differing from that of Noll 13 only by normaliza- 
tion constants: £ 



Zj(r, 0) ■ 



c™R™{r)cos(me) 
c?R?(r)sin(m6) 
c° n R° n (r) 



j even, m ± 0 

j odd, m * 0 » (14) 

m = 0 



where 0 r < 1 and 0 < 6 < 2tt, 

C « [(2 - 5 m0 )(n + DM 1 ' 



(15) 



are normalization constants (the factor fir in the denomi- 
nator is the only difference from the notation of Noll 13 ), 



(n~mV2 



(-l)'(n-s)! 



s! [(n + m)/2 - 5]! [(n - m)/2 - s]\ 



(16) 



are the radial Zernike polynomials, and n and m must 
be positive integers satisfying m ^ n, n - m even. The 
index m is the azimuthal frequency of a given Zernike 
polynomial, and n is its radial degree. The number j is 
a convenient mode-ordering index; it may be verified that 
in the ordering of Noll 13 each valid pair of indices (m, n), 
m ^ 0, corresponds to the pair of consecutive integer 
numbers (j, j + 1), where 



. n(n + 1) ^ 
J = j{m, n) = + m, 



(17) 



whereas each pair (0, n) corresponds to only one num- 
ber j + 1, j = j(0 t n) from Eq. (17) (see Table 1). Impor- 
tantly, this ordering has no gaps; i.e., if we consider the 
set of all Zernike polynomials with radial degree not ex- 
ceeding some integer N, their j indices will constitute the 
set Jn = {1, 2, . . . , js] of all consecutive integers from 1 
to jx, where 

Jn = j(N, N) + 1 = (N + l)(N + 2)/2. (18) 

The choice (15) of normalization constants makes the sys- 
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tern of Zernike polynomials orthonormal with respect to 
the standard scalar product: 



<Z, 



, Zj) « f f Zi(r t d)Zj(r, $)r drdd = 6 U . 
Jo ./o 



(19) 



Let us introduce the vector space Z N of all linear com- 
binations of Zernike polynomials with radial degrees not 
exceeding an integer N: 



X QjZjj Q-j real numbers 



(20) 



(see the example for N = 5 in Table 1). It follows from 
Eqs. (18) -(20) that the dimension dim(Z^) of the space 
Z N is equal to j N : dim(Z^) = (N + l)(N + 2)/2. 

We will also need the subs pace DZ N spanned by all 
diagonal polynomials from Z N : 



VZ N = [ZajZj : Zj GZ N ,m = n] 



(21) 



(see Table 1). Obviously, dim(DZ„) = 2N + 1, as there 
are two diagonal Zernike. polynomials c"r n sin(rc0) and 
c£r n cos(n$) for each radial^legree n * 0 and one diagonal 
polynomial of zero degree, Zo = Cq = ir~ m . Note the 
important equality 



dim(Z N ) - dim(DZ^) « dim(Z^_ 2 ), 



(22) 



which follows directly from the formulas for the dimen- 
sions presented above: (N + 1){N + 2)/2 - (27V + 1) = 
(N - l)N/2. 

Let us consider the restriction (-A)# of the Laplacian 
to the finite-dimensional space Z N . One can easily verify 
that the polynomial (-A)Z; has the same form as that of 
Zj [see Eqs. (14)- (16)], with different radial components 
R?(r) in place of R?(r): 



(n-m)/2 
s-0 



(23) 



We denote by Ker[(-A)jv] the kernel of the operator 
(-A)n, i.e., the set of all functions <f> G Z# mapped by 
(-A)^ into zero. We state that the kernel of (-A)# co- 
incides with the space spanned by the diagonal Zernike 
polynomials: 



Ker[(-A)„] = DZ*. 



(24) 



In other words, the Laplacian of a linear combination of 
Zernike polynomials is equal to zero if and only if this lin- 
ear combination contains only the diagonal polynomials. 
It is easy to see from Eq. (23) that AZ, = 0 for any diag- 
onal polynomial, which implies that Ker[(-A)jv] D DZ;v. 
The opposite inclusion, i.e., the fact that any function <p 
from Z//, such that A<£ = 0, can be expressed as a lin- 
ear combination of diagonal Zernike polynomials, is less 
obvious and is proved in Appendix A. 

Now consider the image space Im[(-A)/v] of the opera- 
tor (-A)^, i.e., the set into which Zs is mapped by the 
Laplacian: Im[(-A)#] = (-A)Z^. We state that the im- 
age of (-A)at coincides with the space spanned by the 
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Zernike polynomials with radial degrees not exceeding 
N - 2: 

lm[(-*) N ] = Z N - 2 - (25) 

In other words, a function is equal to the Laplacian of a 
linear combination of Zernike polynomials with radial de- 
grees not exceeding N if and only if it can be expressed 
as a linear combination of Zernike polynomials with ra- 
dial degrees not exceeding N - 2. To prove Eq. (25), we 
note, first, that according to Eq. (23) the Laplacian de- 
creases the radial degree of any Zernike polynomial by 2. 
It indicates that Im[(-A)/v] C Z#- 2 . A rigorous proof of 
this inclusion is given in Appendix B. Second, the di- 
mension of the image space of any linear operator is al- 
ways equal to the difference between the dimension of the 
whole space, where it is defined, and the dimension of its 
kernel. In the case of (-A)jv we have 

dim{Im[(-A) N ]} - dim(Z^) - dim{Ker[(-A) N ]} 

= dim(Z w ) - dim(DZjv) = dim(Z w . 2 ), 

(26) 

where we used formulas (22) and (24). Thus we see that 
the image of (-A)# is a subspace of Z^- 2) and its dimen- 
sion is equal to the dimension of Z^_ 2 . Therefore the 
vector subspace Im[(-A)#] is equal to the whole vector 
space Z/v-2, and Eq. (25) is proved. 

Let us denote by UZw the subspace of Z^ spanned 
by all nondiagonal polynomials (see Table 1). Because of 
the orthogonality relations (19), any function <f> from Z^ 
can be uniquely represented as a sum of two orthogonal 
components: 

<f> = <f>™ + <£ (1 >, <£ (0) G DZw, <p {1) e UZ* - (27) 

In other words, we have the decomposition of the 
space Z N into an orthogonal sum of two subspaces: 
2tN — DZtf ® UZtf, where © denotes the orthogonal sum. 
Formulas (24) and (25) give us the corresponding decom- 
position of the Laplacian: 

<-«»*»-[; -a](™i)-( z °,)' 

""-""tal-C)" 28 ' 

i.e., operator (- maps the diagonal subspace DZ# into 
zero, and it maps UZ/v one to one (bijective) on Z^-2. 
Hence the restriction of the Laplacian to the subspace of 
nondiagonal Zernike polynomials is invertible, and corre- 
sponding components of the phase can be uniquely found 
from the wave-front curvature. On the other hand, the 
diagonal component 4>f^ cannot be obtained from the 
Poisson equation, and we will need to use the bound- 
ary conditions (9) for its determination. The actual al- 
gorithm is described in Subsection 3.C. 

C. Zernike Decomposition of the Neumann 
Problem [Eqs. (8) and (9)] and Its Solution 

In this subsection we will use the representation (28) of 
the, Laplacian and the Neumann boundary condition (9) to 
derive an algorithm for the unique reconstruction of the 



phase <f> from the wave-front curvature f and the boundary 
slope ip. 

Suppose that the wave-front curvature f obtained 
from intensity measurements is approximated by a fi- 
nite Zernike sum f {N >\ G Z^' with some integer TV' (we 
have explained in Section 2 that the phase solution is 
stable with respect to small errors in f that may occur as 
a result of this approximation). We will be looking for a 
phase <f> such that -A<£ = fa*). -Let AT = N' + 2. Then 
f(N') = f{N-2) £ Zw_2, and, according to Eq. (28), we must 
look for the corresponding phase solution in the space 
Ztf, i.e., find the coefficients 4>j in the representation 
<t> {N) (Rp, 0) - X; e</iV <t>jZj{p, 0), where p = r/R. 

We start with the decomposition of the function <£ (lV ) G 
Zjv in accordance with Eq. (28): 

<f>(N) = <f>(N) + 0(AT), 4>W) = X <f>jZj, 

<C> = I 4>jZ jt (29) 

where the set DJs contains all indices j from J$ - 
{1, 2, .. ., j N ] corresponding to the diagonal polynomials 
Zj and the complementary subset UJ N = J N \DJ N con- 
tains all indices corresponding to nondiagonal ones. Ac- 
cording to the results of Subsection 3.B, we should be 
able to retrieve the nondiagonal component from the 
wave-front curvature fw-2) using the Poisson equation. 
Let us decompose f(N-2) into the Zernike terms, 

f(N-2)(Rp,e)= I fiZi(p t 6) t (30) 
and use Eq. (28) to obtain 

X fiZi = f ( N-2) = {-&)N<f> { (N) = (-A)* I <t>jZj 

= «" 2 I I *>AyZ,-. (31) 

Here Aij = <-AZj, Z ( -> for all i G Jn-i, j £ UJ N> with 
the scalar product ■) defined in Eq. (19). Note that, 
on account of the equality (22) and the definition of the 
set UJtfy the matrix [A ,-_,•] is square with both dimensions 
equal to j^-2 — N(N - l)/2. It is also invertible (non- 
singular) because of Eq. (28). As the Zernike polynomi- 
als Zj are linearly independent, all coefficients at Z, in 
Eq. (31) with the same indices must be equal, which gives 
us the following system of linear algebraic equations: 

Z Aij<t>j = R 2 f h i<=J N - 2 . (32) 

Solving Eq. (32) for we obtain 

<f>j = R 2 I AJiYi, JEUJ N , (33) 

where [AJ; 1 ] = [A,y]~ l is the inverse matrix [it represents 
the operator (-A)^ 1 : Zjv- 2 — • UZ^]. 

Thus we have retrieved the nondiagonal component 
of the phase using Eq. (33), and it remains for us 
to find the diagonal component We will use the 

boundary condition (9) for this purpose. As <f> W) belongs 
to the space Z^, it contains circular harmonics of the or- 
ders m < N. Consequently, its normal derivative at the 
boundary, d r </>(N){R, &), belongs to the space F# spanned 
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by all circular harmonics with azimuthal frequencies not 
exceeding N: 



Fv- 



V(0): VW = vo + I [v' m sin(m0) 



+ ^ cos(m0)] 



(34) 



We assume that the boundary wave-front slope tp obtained 
from the intensity measurements is approximated by a 
finite Fourier sum tf/ iN) E F^: 

N 

= ^0 + I sin(m0) + ^ cos(m0)]. (35) 
Then we can write the boundary condition (9) as 



(36) 



Obviously, d r <f> m (R t $) = a r <0* > 6) + 3 r 0), and 
we can calculate the derivative drfi^R, $), as the <£ ( ^ 
component is already found. Let us introduce the new 
function tj/ iN) (6) = (A<*)(0) - dr^Sotf*, 8) and decompose 
it into the circular harmonics: 



tj>(N)(0) = 4fo + X sin(m0) + ^ cos(m<9)]. 

m-l 



(37) 



On the other hand, by its definition <f>f^ is a linear com- 
- bination of diagonal Zernike polynomials: 

4>w(Rp> e ) = c l<f>o + X Cp" 1 ^ sin(ro0). 

m-l 

+ 0 m cos(m0)]; (38) 
hence its radial derivative can be written as 

d r <f>Z(R, 6) = /T 1 I cZm[<t>' m sm(m$) + 0£ cos(mfl)]. 



(39) 



Finally, as d r <f>l%(R, d) - we obtain from Eqs. (37) 

and (39) the following equations for the diagonal coeffi- 
cients of the phase: 



RK 



cannot be valid [i.e., Eq. (37) cannot be equal to Eq. (39)1, 
unless tp 0 = 0 in Eq. (37). By the conventional formula 
for the Fourier coefficients we have the equality 



1 f 2 " 

= 2^/0 1* 



'imW) ~ dr«,V,(/*. 0)]d0 = 0. (41) 



Applying Green's theorem to Eq. (41), we obtain 

f f f (N - 2) (Rp, 0)pdpde. (42) 

Jo Jo Jo 

If we substitute the Fourier expansion (35) for and 
the Zernike expansion (30) for f {N - 2 ) into Eq. (42), we see 
that the integrals of all the components, except the zero- 
order ones, are equal to zero because of the orthogonality 
of the system of circular harmonics and Zernike polyno- 
mials. Thus condition (42) is equivalent to 



00 = 



2^F 



fo. 



(43) 



where ^0 is the zero-order Fourier coefficient of and f 0 
is the zero-order Zernike coefficient of/". We can also ob- 
tain Eq. (43) from the eneigy-conservation law, substitut- 
ing the Fourier expansion for ^ and the Zernike expansion 
for f in Eq. (10). Therefore Eq. (43) and hence Eq. (41) 
are always true. 

Thus our reconstruction algorithm allows the unique 
retrieval by means of formulas (33) and (40) of all phase 
aberrations with radial degree not exceeding N, except 
the piston mode. In Subsection 3.D we will derive an 
explicit matrix formula corresponding to this algorithm, 
which expresses the reconstructed phase as a function 
of Zernike coefficients of the wave-front curvature f and 
Fourier coefficients of the boundary slope tp. 

D. Phase Retrieval Matrix 

We consider the operator A N defined by the Laplacian and 
the Neumann boundary condition in the space Z Ni A N = 
{(-A)w, B N ), with {B<f>){6) = d r <t>{R, 0), and write the ana- 
log of decomposition (28) for it: 



(44) 



m = 1, 2, .... N . (40) 



Formula (40) gives us the values of Zernike coefficients 
of the diagonal phase component <f>l%, except the zero 
coefficient <f> 0 . Applying formula (17) with m = n, we 
derive that each integer m in Eq. (40) corresponds to a 
pair ofy indices (j 9 j + 1) = [m(m + 3)/2, m(m + 3)/2 + 1] 
of a pair of diagonal polynomials. Hence, when m runs 
through the sequence 1, 2, . . . , N, these pairs (j, j + 1) 
exhaust the whole set DJn, except the zero-order Zernike 
polynomial [note that m starts from unity in Eqs. (39) 
and (40)]. 

The last exception has two important consequences. 
First, the impossibility of determining the coefficient <f> 0 
means that we can reconstruct the phase only up to an 
arbitrary additive constant. This, however, is a well- 
known general property of the Neumann problem already 
discussed in Section 2. Second, boundary conditions (36) 



where B^ 




(45) 



DZ N - F N and B ( „,: UZ N — are the 
restrictions of the boundary operator {B<f>){6) = S r <f>(R, 9) 
to the corresponding subspaces. Let us also introduce 
the spaces t iN) and DZ,^, obtained from F (A0 and DZ )/V| , 
respectively, by the removal of the zero-order (constant) 
component and define the corresponding operators: : 
DZ# — ► ¥ N and B [N) : \JZ N ^F N . We must remove the 
zero-order component to make the operator A N invertible 
according to Subsection 3.C. Let us rewrite Eq. (45) in 
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the Zernike and Fourier spaces without constant compo- 
nents: 



7 [ 4>iN) 



fi(O) 



fid) 
tftf 



0 



L <f>W) I 



, p<i) 



= (/"')■ 

\/(N-2)J 



(46) 



Now one can easily find the inverse operator A^ 1 by ma- 
trix algebra: 



An ~l 0 
A ^(z?- 2 ) = (uzl)' 



(47) 



VMM- 2)/ 



[A%] ^(N) 



" [^Ar)]" 1 ^)!"^^)] V(tf-2> 

[~^(A/)]~ I /'(iV-2) 



(48) 



Formula (48) is an operator representation of the recon- 
struction algorithm described in Subsection 3.C. Each 
term in Eq. (48) has already been defined. Operator 
[-A (A0 ] _1 is defined in Eq. (33). It is represented by 
the matrix R 2 A]* applied to the vector of Zernike coeffi- 
cients of /W-2). Operator [B^]* 1 is defined in Eq. (40). 
Its action on a function tp {N) from F iN) is equivalent 
to the multiplication of its Fourier coefficients $' m9 
by the constants R/mc™. Operator B$> is the restric- 
tion of the boundary operator {B </>)($) = d r <f>(R, 0), &$ } : 
UZ;v — FV The inversion (48) admits a generalization 
onto the infinite-dimensional case, 24 but we will not need 
it in this paper. 

Thus we have obtained the explicit formula (48) for the 
j operator A^ 1 , which maps each pair {^(^ )f fa -2)} of the 
r wave-front curvature fa -2) £ Z^_ 2 and boundary slope 
ify{N) £ F^ onto the unique (up to an additive constant) 
phase solution 4>w) — <f>(N) + <f>(N) S Z/y such that 
-&<f>W) = fa-2) and d r <f> iN) (R,6) = tj; {m (6). Formula 
(48) is the central result of the present paper. 

Now we will give the analytical form for the ma- 
trix formulas (46)-(48) convenient for applications and 
present a numerical example for N = 4. Let us start 
from the matrix A#: 



- fA u A 12 ~ 



R-'B^ R- l B?j 

0 R~ 2 A; 



(49) 



where 



(H 



c = (5D 



i odd 



where i = 1, 2, . . . , 2N, j G DJ S \{1] for B\j\ and j G UJ N 
for Bij\ Using Eq. (50) and the definition (14) -(16) of 
the Zernike polynomials, we obtain 



R- l B?} } -R 



r-'b'^r- 1 



i = 1, 2, .... 2N, jsDJ N \{l} 9 m=m(j), (52) 

vjSik, i = l,2 2^)6%, 

(53) 

2m(j) - 1 j odd 
2m( i /') 7' even 



(54) 



where <r,. = ZiV^ yl. m (» " 2 *) [see Eqs. (15) and (16) 
for the definition of the Coefficients and y* m ], 8 ik is the 
Kronecker symbol, and m(j) is the azimuthal frequency 
of Zj. The block A 22 « R " 2 Ay is the (AT - l)N/2 x (N - 
l)N/2 square matrix with elements 

R- 2 Aij-(-*Z j9 Zi) 

= [ [ (~A)Zj(r/R s 6)Zi(r/R t 6)rdrdd, (55) 
Jo Jo 

where i G J N - 2 and j G UJ N . The following formula is 
derived from the definition (14)-(16) of Zernike polyno- 
mials and formula (23) for the Laplacian of Zernike poly- 
nomials: 

(n'-ml/2 (n-m)J2 

Aij = S mm .2[(n + 1Kb' + l)] m X I TikriU 

s'-0 s-0 

m 2 - (n - 2s) 2 



n f + n - 2s' - 2s 



(56) 



where i G J N ~ 2 , j G £A/*, J = j(m, n), and ( = j{m\ n'), 
as in Eq. (17), and the coefficients y* ttn are defined in 
Eq. (16). Note that although formulas (52), (53), and (56) 
give explicit analytical expressions for the elements of 
the matrix A ( jv), it is sometimes more convenient to use 
formulas (50) and (55) for the practical calculations. In 
Table 2 we present an example of matrix A ( jvj for N — 4. 
The inverse matrix A ( 7J, consists of the blocks 



a n - 



An 1 



0 



-RKB%]-*B$A? 
R 2 A 



(57) 



In Eq. (49) A n = R- x B{j is a square 2N X 2N matrix, 
and A 12 = J R- 1 B t -j , is a rectangular 2N x (N - l)N/2 
matrix with elements 

f 2n 

R-'Bij - (BZ j9 Fi) r - R' 1 (d p Zj)(l 9 0)F i (0)d6 , (50) 
Jo 



In Eq. (57) A\} is a square 2N X 2N matrix, All is 
a rectangular 2N x (N - l)N/2 matrix, and A 22 is an 
{N - l)N/2 x (N - l)N/2 square matrix. In order to 
obtain the elements of A^ )t one needs to find only the in- 
verse matrices A"; 1 = (A 0 )* 1 and [B\f]~ l . The elements 
of these inverse matrices can be easily calculated: 
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Table 2. Matrix A (4) for R - 1° 



7 





2 


3 


5 


6 


g 


10 


14 


15 


4 


7 


8 


1 1 


12 


13 


1 


0 


2 


0 


0 


0 


0 


0 


0 


0 


14>/2 


0 


0 


0 


0 


2 


2 


0 


0 


0 


0 


0 


0 


0 


0 


0 


14>/2 


0 


0 


0 


3 


0 


0 


2>/6 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


10VT0 


4 


0 


0 


0 


2>/6 


0 


0 


0 


0 


0 


0 


0 


0 


loyio 


0 


5 


0 


0 


0 


0 


6y/2 


0 


0 


0 


0 


0 


0 


0 


0 


0 


6 


0 


0 


0 


0 


0 




0 


0 


0 


0 


0 


0 


0 


0 


7 


0 


0 


0 


0 


0 


0 


0 


4VT0 


0 


0 


0 


0 


0 


0 


8 


0 


0 


0 


0 


0 


0 


4>/U> 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


0 


0 


-873 


0 


0 


-24V5 


0 


0 


2 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


-24^2 


0 


0 


0 


3 


0 


0 


0 


0 


0 


0 


0 


0 


0 


-24v/2 


0 


0 


0 


0 


4 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


-16>/l5 


0 


0 


5 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


-16VT5 


6 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


-16^ 


0 



°The rules separate the blocks A M . The indexing system is described in the text [formulas (49) and (52>-(56)]. 



Table 3. Phase Retrieval Matrix A ( ~ 4 | for R = l a 

i 



j 1 


2 


3 


4 


5 


6 


7 


8 


1 


2 


3 


-* 4 


5 


6 


2 0 


1/2 


0 


0 


0 


0 


0 


0 


0 


7/24 


0 


0 


0 


0 


31/2 


0 


0 


0 


0 


0 


0 


0 


0 


0 


7/24 


0 


0 


0 


5 0 


0 


1/(2 J6) 


0 


0 


0 


0 


0 


0 


0 


0 


0 


5/48 


0 


6 0 


0 


0 


1/(2 V6) 


0 


0 


0 


0 


0 


0 


0 


0 


0 


5/48 


9 0 


0 


0 


0 


l/(6>/2) 


0 


0 


0 


0 


0 


0 


0 


0 


0 


10 0 


0 


0 


0 


0 


1/(6 V2) 


0 


0 


0 


0 


0 


0 


0 


0 


14 0 


0 


0 


0 


0 


0 


0 


1/(4^) 


0 


0 


0 


0 


0 


0 


15 0 


0 


0 


0 


0 


0 


l/(4>/l0) 


0 


0 


0 


0 


0 


0 


0 


4 0 


0 


0 


0 


0 


0 


0 


0 


-1/(8V3) 


0 


0 


1/16 


0 


0 


7 0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


-1/(24^) 


0 


0 


0 


8 0 


0 


0 


0 


0 


0 


0 


0 


0 


-1/(2W2) 


0 


0 


0 


0 


11 0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


-1/(16715) 


0 


0 


12 0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


-1/(16715) 


13 0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


-l/( 16^15) 


0 



°The rules separate the blocks A^J. The first 14 nonconstant Zernike aberrations (j = 2, 3, .. . , 15) of a phase can be retrieved by the application 
of this matrix to the vector {iPn), f<2)) of the Fourier coefficients of the boundary slope and' the Zernike coefficients of the curvature of the wave front. 
The indexing system is described in the text [formulas (57)-(60)]. 



£ = 1, 2 2N, j e DJ N \[l} t m = m{j) , (58) 

(A^)ji = R 2 (-ir j Mij/det A, i G J N . 2 , j e UJ N , 

(59) 

where det A is the determinant of the W - l)N/2 x {N - 
l)N/2 matrix A;; with the elements defined in Eq. (56) 
and Mij is the determinant of the matrix obtained from 
Aij by removing the ith row and the jib. column. Thus 
we defined all the elements of the diagonal blocks A7i 
and A^ 1 - The elements of the nondiagonal block All 
can now be obtained according to the rules of matrix 
multiplication: 

(At})* = -R 2 1 KAllhkBu'iA^hi, 
k 1 

j e DJ N \{i) t ieJ N -2- (60) 

Formulas (58)- (60) give explicit analytical expressions 
for the elements of the matrix Aj/ ; however, in practice, 



it may be easier to obtain the whole matrix A^Jj by the t .} 
matrix inversion of A ( ao. 

The phase £ Z ( ,V) is retrieved by the direct mul- 
tiplication of this matrix A^ 1 and the vector g {N ^ = 
(^1, f{N-2)} of the Fourier coefficients of boundary slopes 
tfrw) E and the Zernike coefficients of the wave-front 
curvature f{N-2\ £ Z/v_ 2 : 




(61) 



As can be seen from Eqs. (57)-(60), the phase retrieval 
matrix A(~^, does not depend on experimental data and 
can be calculated for a given W once and for all with arbi- 
trary precision. In fact, using algebraic, rather than nu- 
merical, calculations with the help of such mathematical 
software as, for example, MATHEMATICA, one can calculate 
the matrix A,7Jj with absolute precision. We performed 
this algebraic calculation for N — 4, with the result pre- 
sented in Table 3. 
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It is worth noting that the matrix in Table 3 is very 
sparse and that the values, of its nonzero elements are 
decreasing as the matrix indices increase. This indi- 
cates that the retrieval of phase with this matrix should 
be quite stable. We still expect larger errors in the re- 
trieval of the diagonal (with m = n) Zernike aberrations 
compared with those arising from the nondiagonal terms 
(with m # n) t because of the necessity to distinguish 
the contribution of the wave-front boundary slopes from 
that of the wave-front curvature. As can be seen from 
Eq. (61), if the uncertainties in & tNi are greater than the 
uncertainties in fa-2i, then the uncertainty in the recon- 
struction of the diagonal phase component <j>l% will also 
be greater than that in the nondiagonal component </>lll ir 
as the latter depends only on the wave-front curvature 
fw-2). A full discussion of the stability phase retrieval 
with this approach is beyond the scope of the present 
paper, and a thorough study of stability is currently 
under way. 

4. SUMMARY 

In this paper we have examined the solution of the 
transport-of-intensity equation (TIE) from the perspective 
of a Zernike polynomial decomposition. We have found 
a direct matrix relationship between a Zernike polyno- 
mial decomposition of the intensity derivative and the 
Zernike decomposition of the phase. This relationship 
gives a precise description of the influence of individ- 
ual Zernike aberrations of the phase on the evolution 
of the intensity distribution in the wave front near a 
uniformly illuminated circular aperture, making phase 
retrieval extremely straightforward. The derived opera- 
tor solution (48) to the Neumann boundary problem for 
the TIE and its matrix representation (61) give a unique 
phase 0^ E Z ( ^> for any given wave-front curvature in- 
side the aperture fa -2) E %w-2) and the boundary slope 
*lt\N\ E ¥ {N) . In practice, <j> {N) is obtained as a sum of 
two orthogonal components and 0$,, the first 

containing all diagonal Zernike aberrations of the phase 
and the second containing all the nondiagonal ones. The 
component <£ (iV , depends only on the wave-front curva- 
ture inside the aperture. The diagonal component 
depends on the wave-front slope at the boundary as well 
as on the boundary values of . We have also derived 
explicit analytical expressions [Eqs. (57)-(60)] for the 
phase retrieval matrix A^ 1 corresponding to the operator 
solution (48) and presented a numerical example of such 
a matrix for N = 4. The phase <f> {Ni E Z (N) can be re- 
trieved by the direct multiplication of the matrix A^ with 
the vector g (N) = fa- 2 )\ of the Fourier coefficients of 
boundary slopes and the Zernike coefficients of the wave- 
front curvature. The matrix A^ 1 does not depend on 
the experimental data and so needs to be calculated only 
once for a given N. Thus there is the prospect of a very 
rapid and precise phase retrieval algorithm with this ap- 
proach. We are planning to present results of computer 
simulations with this algorithm in a subsequent paper. 

APPENDIX A 

Here we prove that the kernel of the Laplacian in the. 
Zernike space Z lV coincides with the subspace DZjv 



spanned by the diagonal Zernike polynomials, i.e., 
Ker[(- A);v] = DZ jV . It was shown in Subsection 3.B that 
&Zj = 0 for any diagonal polynomial, so Ker[(-A)^] D 
DZ;v. Hence it is sufficient to prove the opposite inclu- 
sion, 

Ker[<-A)*]C DZ jV , (Al) 

i.e., that any function £ from Zy, such that A<£ = 0, can 
be expressed as a linear combination of diagonal Zernike 
polynomials. Let us consider an arbitrary function <f> 
from Z N : 

<t> = X <t>jZj = X Zr a [6' nm sin(mtf) + 4>'n„ cosfm*)], 

(A2) 

where the first sum is over n < AT and the second is over 
such m that n - m is positive and even. Then 

-A* - I Km 2 - n 2 )r«- 2 [<l>' nm sin{m6) + <j>L cos(mtf)]. 

n m 

(A3) 

All the monomials r n ~ 2 sin(m0) and r n ~ 2 cos(m#) in 
Eq. (A3) are linearly independent, so A<£ = 0 if and only 
if all the coefficients in*Eq. (A3) are equal to zero: 

(m 2 - n 2 )4>' nm = 0, (m 2 - n 2 )tf M = 0 . (A4) 

But Eqs. (A4) imply that all coefficients <t>' nm and <j>$ m with 
n * m must be equal to zero, so in the expansion (A2), only 
diagonal coefficients <b' mm and ^ m may be nonzero, i.e., 
the function <f> belongs to the diagonal subspace DZ^ and 
relation (Al) is proved. 

APPENDIX B 

Here we prove that the Laplacian maps the Zernike space 
Z# into the space Z/v- 2 , i.e., 

Im[(-A)iv]CZ N -2. (Bl) 

In other words, we must prove that for any function xf> 
from Ztf its Laplacian A<£ can be represented as a lin- 
ear combination of some Zernike polynomials with radial 
degrees not exceeding N - 2. As Arf» can be expressed 
by formula (A3), it is sufficient to prove that each mono- 
mial r n ~ 2 sin(m0) and r n ~ 2 cos(m0) in Eq. (A3) belongs to 
Zn-2- Importantly, m < {n - 2) for all terms in Eq. (A3), 
because the terms with m — n from Eq. (A2) are mapped 
by the Laplacian into zero. We consider only monomials 
with sin (the proof for the monomials with cos is iden- 
tical). Thus we are going to prove that any monomial 
r l sin(m0), with / < n - 2, / - m nonnegative and even, 
belongs to Z^_ 2 * It is then sufficient to prove that for any 
m < N - 2 the radial monomial r', / = m + 2k, I < AT - 2, 
can be represented as a linear combination 

r l = X <T*R?ir) (B2) 

/IS/ 

of radial Zernike polynomials with fixed m. We will 
prove it by induction by k. For k = 0 we have / = m 
and Eq. (B2) is obvious, as r m = R% [see Eq. (16)]. Let 
us assume that we have proved Eq. (B2) for / = m + 
2(k - 1) and then prove it for / = m + 2k. Using the 
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definition (16) of Zernike polynomials, we can write 



Ym*ik.,„r m 



I 



But by the induction assumption all monomials under the 
summation sign in Eq. (B3) can be represented in the 
form (B2) [as their indices do not exceed m + 2{k - 1)], 
and R'n t +2k obviously has the form ( B2) with m + 2k — I < 
N - 2. Therefore r" ,+2 * can be represented as the linear 
combination (B2), and hence Eq. (Bl) is proved. 
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